Quantum path integral simulation of isotope effects in the melting 

temperature of ice Ih 
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The isotope efTect in the melting temperature of ice Ih has been studied by free energy calcula- 
tions within the path integral formulation of statistical mechanics. Free energy differences between 
isotopes are related to the dependence of their kinetic energy on the isotope mass. The water sim- 
ulations were performed by using the q-TIP4P/F model, a point charge empirical potential that 
includes molecular flexibility and anharmonicity in the OH stretch of the water molecule. The 
reported melting temperature at ambient pressure of this model (T=251 K) increases by 6.5±0.5 
K and 8.2±0.5 K upon isotopic substitution of hydrogen by deuterium and tritium, respectively. 
These temperature shifts are larger than the experimental ones (3.8 K and 4.5 K, respectively). In 
the classical limit, the melting temperature is nearly the same as that for tritiated ice. This unex- 
pected behavior is rationalized by the coupling between intermolecular interactions and molecular 
flexibility. This coupling makes the kinetic energy of the OH stretching modes larger in the liquid 
than in the solid phase. However the opposite behavior is found for intramolecular modes, which 
display larger kinetic energy in ice than in liquid water. 

PACS numbers: 64.70.dj, 82.20.Wt, 64.70.D-, 65.20.De 
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I. INTRODUCTION 

Quantum mechanical effects associated to the nuclear 
mass play a significant role in the properties of liquid 
water and ice. Experimental evidence is provided by the 
isotope dependence of the equilibrium properties of wa- 
ter. At ambient pressure the melting point at 273.15 
K increases by 3.8 K and 4.5 K after isotopic substitu- 
tion of hydrogen by deuterium or tritium, respectively. 
An even larger isotope effect is found in the tempera- 
ture (r=277.13 K ) of maximum density (TMD), that 
increases at ambient pressure by 7.2 K in heavy water 
and by 9.4 K in tritiated water. Such behavior can not 
be described by classical statistical mechanics, as in this 
limit the atomic masses do not affect the phase diagram 
of a substance or the equilibrium structure of a liquid. 
The importance of quantum effects related to the atomic 
masses in water might be expected by the presence of the 
lightest atom. However, the larger oxygen mass is also 
the origin of significant quantum properties of even prac- 
tical relevance. For example, the isotopic composition of 
the annual layers of ice accumulated in the Antarctica 
has provided an indirect measure of the temperature of 
our planet over the last 400,000 years. The reason is that 
the vapor pressures of H2^^0, H2^*0, and D2^^0 are dif- 
ferent, and then the isotopic composition of ice results to 
be a function of the temperature at which it precipitated. 
Thus, isotope analysis in ice provides an historical record 
for climate change in the pastii 

Computer simulation of water in clusters and con- 
densed phases has attracted a lot of interest since the 
pioneering work of Barker and Watts^ and Rahman and 
Stillinger'^ using rigid nonpolarizable models for the wa- 
ter molecule. Since then a lot of effort has been invested 
in the development and refinement of empirical potentials 
for both water and ice simulations. In fact, there appears 



an embarrassing variety of them in the computer simula- 
tion literature. The most employed models assume a rigid 
geometry of the water molecule, some include molecu- 
lar flexibility either with harmonic or anharmonic OH 
stretches, and another group deals explicitly with po- 
larizability effects.'' Moreover, in some cases slight mod- 
ifications of the potential parameters are proposed for 
their use in quantum simulations, to avoid overcounting 
of quantum effects if the model parameters were first fit- 
ted against experimental data by classical simulations. 
Besides, there is an increasing number of water simula- 
tions using ab initio density functional theory (DFT)iii^ 
However, the H-bond network, with a strength between 
weak covalent and van der Waals interactions, seems diffi- 
cult to be described with presently available energy func- 
tional. As a result, some properties of ice may be poorly 
reproduced by DFT simulations, e.g., its melting temper- 
ature can be overestimated by more that 130 KM- 

The melting point of ice at ambient pressure has been 
determined for the most common rigid models within the 
classical limitji^ It was found that the T1P4P model, with 
melting point at T = 232 K, results superior to other 
models in the sense that correctly predicts that ice Ih is 
the stable phase at ambient pressure, while the predic- 
tion of the other rigid models (SPC, SPC/E, T1P3P and 
T1P5P) was ice 11. An improved parametrization of the 
rigid model (T1P4P/2005) displays a melting point at 251 
Kiii Quantum simulations of phase coexistence are less 
common than their classical counterparts. An exception 
is the work of Habershon et al.,^^ who have developed a 
fiexible water model (q-TlP4P/F) by adding to the rigid 
T1P4P/2005 potential the intramolecular flexibility with 
the help of OH Morse-type stretches. The model was 
parametrized on the basis of quantum path integral (PI) 
simulations. Its melting point at 251 K and ambient pres- 
sure, was derived by direct coexistence PI simulations of 
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the water-ice interface. Moreover, in the classical limit 
the melting point was found just 8 K above the quantum 
result. This temperature shift was considered consistent 
with the experimental difference of 4 K between the melt- 
ing points of H2O and D20.^^ For the point charge flex- 
ible q-SPC/Fw modeli^ the classical melting point was 
however found about 27 K higher than the quantum re- 
sult of 195 K.^^ This difference in the quantum correction 
of both models might be originated from the description 
of the intramolecular OH stretches, i.e., anharmonic (q- 
TIP4P/F) versus harmonic (q-SPC/Fw) OH vibrations. 

Although the determination of quantum corrections 
to classical melting points is interesting because it al- 
lows us to quantify the systematic error of treating water 
molecules as classical entities, they do not represent any 
kind of measurable property. There is no way to perform 
measurements of the phase behavior of water in the clas- 
sical limit. In this respect, the calculation of the isotope 
effect in the melting point of ice offers the advantage of 
being directly comparable to experimental data, provid- 
ing a better test of the capability of the water model. We 
are not aware of previous computer simulations of this 
isotope effect. However, by the quantum cluster equi- 
librium theory, that calculates equilibrium properties by 
extending standard quantum statistical thermodynamics 
of chemical equilibrium to the analogous equilibrium be- 
tween molecular clusters, it was estimated that the melt- 
ing point of D2O is shifted by about 2 K towards higher 
temperatures with respect to light water. Isotope ef- 
fects have been studied by PI simulations in many other 
equilibrium properties of water. Kuharski and Rossky 
found that the liquid H2O is less structured than D20.^- 
The explanation was formulated in terms that the quan- 
tum effect associated to the lower isotope mass results 
in a less structured H-bond network and a less tightly 
bound liquid. ^^'^^ Other simulations of isotope effects fo- 
cused on the TMDiii^, the diffusion coefficient the 
heat capacity ,^^'^^ and the infrared spectra^° of water. 

In this paper we present a PI simulation of the isotope 
effect in the melting temperature of ice Ih at ambient 
pressure. The q-TIP4P/F model has been chosen for the 
simulations because it is an anharmonic flexible potential 
whose normal melting point has been already established 
by quantum PI simulations;^- Thus, assuming the equal- 
ity of the Gibbs free energy, G, of ice Ih and water at the 
melting point, the isotope effect will be calculated from 
the dependence of G with isotope mass and temperature. 
Solid-liquid coexistence will be also studied in the classi- 
cal limit. The calculation of G will be performed using 
adiabatic switching (AS) (Ref. 23) and reversible scal- 
ing (RS) (Ref. [13) approaches, that are based on algo- 
rithms where a Hamiltonian parameter (e.g., an atomic 
mass) or a state variable (e.g., the temperature) changes 
along a non-equilibrium simulation run. The capability 
of both AS and RS methods to calculate free energies in 
the context of PI simulations has been recently analyzed 
in the study of the phase diagram and isotope effects of 
neoni^i2^ 



The structure of this paper is as follows. In Sec. Illlwe 
present the computational conditions employed in the PI 
simulations as well as the techniques used to evaluate 
the free energy as a function of the isotope mass and 
temperature. In Sec. IIIIl selected results of quantum 
simulations are compared to available data of Habershon 
et alr^ and also to the classical limit as a check of the 
employed computational conditions. In particular, quan- 
tum and classical radial distribution functions (RDFs) of 
the liquid phase are presented in Subsec. IIII Al while the 
quantum and classical TMD of water at ambient pres- 
sure is summarized in Subsec. IIIIBI The isotope effect 
in the melting temperature of ice is the focus of Sec. 
IIVI Results obtained for D2O and T2O are compared to 
available experimental data in Subsec. IIV Al The clas- 
sical limit is presented in Subsec. IIVBI The calculated 
isotope effects are rationalized by a discussion of the mass 
dependence found for the kinetic energy (KE) in Subsec. 
IIV CI Finally, we summarize our conclusions in Sec. |Vl 



II. METHODOLOGY 

A. Computational conditions 

In the PI formulation of statistical mechanics the parti- 
tion function is calculated through a discretization of the 
integral representing the density matrix. This discretiza- 
tion defines cyclic paths composed by a finite number L 
of steps, that in the numerical simulation translates into 
the appearance of L replicas (or beads) of each quantum 
particle. Then, the implementation of PI simulations 
relies on an isomorphism between the quantum system 
and a classical one, derived from the former by replac- 
ing each quantum particle (here, atomic nucleus) by a 
ring polymer of L classical particles, connected by har- 
monic springs with a temperature- and mass-dependent 
force constant. Details on this computational method are 
given elsewhere, ^^"'^'^ The configuration space of the clas- 
sical isomorph can be sampled by a molecular dynam- 
ics (MD) algorithm, that has the advantage against a 
Monte Carlo method of being more easily parallelizable, 
an important fact for efficient use of modern computer 
architectures. Effective reversible integrator algorithms 
to perform PI MD simulations in either NVT on NPT 
ensembles {N being the number of particles, V the vol- 
ume, P the pressure, and T the temperature) have been 
described in detail in Refs. Isil - fs^ . Both isotropic and full 
cell fluctuations were programmed for the NPT ensem- 
ble. All calculations were done using originally developed 
software and parallelization was implemented by the MPI 
library. •^^ 

Simulations of water were performed on cubic cells con- 
taining 300 molecules, assuming periodic boundary con- 
ditions. Ice simulations included 288 molecules in an or- 
thorhombic simulation cell with parameters (4a, 3\/3a, 
3c), with (a,c) being the standard hexagonal lattice pa- 
rameters of ice Ih. A proton disordered ice structure 
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where each O atom has two chemically bonded and two 
hydrogen bonded H atoms with nearly zero dipole mo- 
ment of the simulation cell was generated by a Monte 
Carlo procedure. Cell fluctuations in the extended dy- 
namics of the NPT ensemble were isotropic for the liq- 
uid phase and flexible for the solid simulation cell. The 
point charge, flexible q-TIP4P/F model was employed 
for the simulations.^^ The Lennard- Jones interaction be- 
tween oxygen centers was truncated at rc=8.5 A, and 
standard long-range corrections were computed assuming 
that the pair correlation function is unity, g(r) = 1 for 
r > Tc, leading to well-known corrections for the pressure 
and internal energies)^ Long-range electrostatic interac- 
tions and forces were calculated by the Ewald method, 
and the calculation was speeded up by allowing the real 
and reciprocal space sums to be performed in parallel. 
The Gaussian smearing parameter in the Ewald method 
was set to 0.465 A, so that both real and fc-space sums 
require similar computation time. 

To have a nearly constant precision in the PI results 
at different temperatures, the number of beads L was 
set as the integer number closest to fulfill the relation 
LT=6000 K, i.e., at 300 K the number of beads was 
L — 20. The classical limit is easily achieved within the 
PI algorithm by setting L=l. The staging transforma- 
tion for the bead coordinates was employed for the quan- 
tum simulations. Temperature was controlled by chains 
of four Nose-Hoover thermostats coupled to each of the 
staging variables, and in the case of the NPT ensem- 
ble an additional chain of four barostats was coupled to 
the volumei^ To integrate the equations of motion, a re- 
versible reference system propagator algorithm (RESPA) 
was employed^^ For the evolution of thermostats and 
harmonic bead interactions a time step dt = At /A was 
used, where At is the time step associated to the calcula- 
tion of the q-TIP4P/F forces. A value of At=0.3 fs was 
found to provide adequate convergence. The virial esti- 
mator was employed for the calculation of the KE^j^ 
and the pressure estimator was identical to that used 
in a previous work.'*" Typical runs consisted of 5 x 10"* 
MD steps (MDS) for equilibration, followed by runs us- 
ing between 5 x 10'^ and 4 x 10^ MDS for calculation of 
equilibrium properties. The longest runs were required 
for the liquid phase, in particular for the calculation of 
the water density as a function of temperature. 



B. Relative free energy 

The thermodynamic integration is a standard tech- 
nique that allows us to obtain free energy differences by 
the calculation of the reversible work needed to change 
the original system into a reference state of known free 
energy4i"— The Hamiltonian is switched along a path 
that connects both systems and a set of equilibrium sim- 
ulations are performed at several points of this path. 
The AS method is an alternative procedure where the re- 
versible work is estimated by slowly changing the system 



Hamiltonian along a single non-equilibrium simulation 
run.'^'^ The RS algorithm was formulated to obtain free 
energies as a function of a state variable, typically T or 
P. In this slow (adiabatic) change of the state 

variable is performed along the non-equilibrium simu- 
lation run.^* Both AS and RS methods have been re- 
cently applied to study the phase diagram of neon by PI 
simulations. ■^^'■^^ Here we summarize how these methods 
fit into our calculation of the isotope effect in the melting 
temperature of ice. 

The first step is to obtain the free energy of each phase 
(solid and liquid) at some appropriate reference point. 
The melting temperature of the q-TIP4P/F model for 
H2O at normal pressure was estimated by Habershon et 
al. by PI simulations.''^^ Hence this state point {Tji=251 
K, Pr=1 atm) is a reference state where the Gibbs free 
energies of the solid, Gs, and the liquid, G/, are identical. 
However, as the melting point was determined by direct 
coexistence, the actual value of the free energy remains 
unknown. This is not a limitation, as only free energy 
differences with respect to the reference state are needed 
for our purpose of calculating isotope shifts. In particu- 
lar, free energy is required as a function of the hydrogen 
isotope mass, mp, and temperature. Thus, without loss 
of generality, we set an arbitrary zero for the entropy, Sq, 
so that we will calculate at ambient pressure relative free 
energies, Gr, defined as 

GR{T,mF)^GiT,mF)-TSo . (1) 

The pressure dependence of G is omitted here as it is a 
constant (P=l atm) for all simulations in this work. An 
alternative to our choice of setting an arbitrary zero of en- 
tropy would be to set an arbitrary zero for the Gibbs free 
energy. Both choices would be physically equivalent in 
the sense that they would lead to identical results for the 
phase coexistence temperature. However, the advantage 
of setting a zero of entropy is that then the temperature 
dependence of Gr is given by an expression identical to 
that valid for G (see Sec. IIIEp . From the last equation 
it is obvious that, at a given temperature, the coexis- 
tence condition of equal free energies of solid and liquid 
phases, Gs=Gi, implies also that the relative free ener- 
gies are equal, Gr^s=Gr^i. The arbitrary zero of entropy 
is formally defined so that the relative free energy of ice, 
Gr^s, (and liquid water Gr^i) is zero at the reference state 
point, {Tr, Pr), for water molecules made of the isotope 
■"■H with mass mn, 

GR,s{TR,mH) = GRj{TR,mH) = Q . (2) 

Our task now is to find the coexistence condition for the 
other hydrogen isotopes at ambient pressure, i.e., the 
temperature, Tm, that satisfies equality in the free en- 
ergies of both solid and liquid phases, 

GR^s{Tm,mF) = GR,i{Tm,mF) , (3) 

for an isotope mass, m^?, corresponding either to deu- 
terium (m£)) or tritium (m^) atoms. 
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C. Isotope effect on the free energy 

If temperature is held constant (say at the reference 
value T/j), free energy differences as a function of the 
isotope mass are the same, whether calculated with G or 
Gr [see Eq. ([T])]. Thus, for simplicity, our derivation is 
presented here by using G and omitting the indication of 
its T dependence. Assuming that the Gibbs free energy, 
G{mH), of a system made of the isotope ■'■H is known, 
then the unknown free energy, G{mp), of a system ob- 
tained by substituting -"-H by the isotope of mass mp can 
be calculated by the expression 



Gimp) = Gimn) + 



dGjmi) 
dmi 



dmi 



(4) 



We recall the thermodynamic relation between deriva- 
tives of the Gibbs free energy and the Hamiltonian, H, 



dG{mi) 
dmi 



dH 
dmi 



NPT 



where the angle brackets represent an ensemble average. 
The Hamiltonian associated to the atomic nuclei in the 
water phase can be expressed as 



H = K{mo) + K{mi) + V 



(5) 



where V is the potential energy operator. K{mo) rep- 
resents the sum of the KE operators of all the oxygen 
nuclei in the system, while K{mj) is the corresponding 
sum for the hydrogen isotopes. The mass dependence in 
the Hamiltonian operator appears only in the KE oper- 
ator, thus 



dH 
dmi 



mi) 



(K{mi) 



NPT 



NPT 



\ d 



'mi 



NPT 



mi 



(6) 

where the last equality follows from the fact that the 
mass mi appears in the KE operator as a mj^ factor. 



K(mi) ) is the ensemble average of the sum of the 

/ NPT 

KE of all the hydrogen isotopes in the system. This quan- 
tity is readily obtained in PI simulations by the virial es- 
timator of the KE, K{mi). Inserting the last result into 
Eq. (jl]), one gets 



G{mF) = G{mH) 



{K{mi)) 



NPT 



mi 



dmi , (7) 



where {K(mi)) j^pj, is the ensemble average of the virial 
KE estimator for the hydrogen isotope. The last equa- 
tion shows that the change in the KE as a function of 
the isotope mass determines the free energy difference 
between two isotopes. We have implemented this free en- 
ergy evaluation by the AS methodjSS The isotope mass. 



TO/, is changed from the initial value {mn) to the final 
one (mp) in a single non-equilibrium simulation. It is 
convenient to change the integration variable, to/, to the 
dimensionless parameter. A/, 



> mn 

Ai — 

mi 



(8) 



so that the integration limits for A/ are A// = 1 and 
\p — mnjmp. For a simulation run consisting of a 
total number J of MDS, the actual value of the isotope 
mass, TO/, at the /'th simulation step (/ = 1,. . . , J), is 
determined by substituting in Eq.([8]) the following value 
of A/ 



A/ = 1 + (/ - 1)AA , 



(9) 



where A A — (\p ~ \h)/{J —1). The Gibbs free energy 
as a function of the isotope mass to/, is discretized for 
the simulation steps / > 1 as 



G(to/) = G(to/_i) + 



1 (K{mi) K{mi_i) 



A/ 



A 



/-I 



AA 



(10) 

where K{mi) is the virial estimator of the KE of the hy- 
drogen isotope at the simulation step /. The change of 
A/ at each simulation step implies to update the corre- 
sponding isotope mass, to/, as well as the spring constant 
for the harmonic coupling between beads of that isotope. 
A convenient application of the AS method is to perform 
two independent non-equilibrium simulations for a given 
free energy determination, where the initial and final in- 
tegration limits are interchanged, so that the reversible 
path between the initial and final integration points is 
run in both directions. The free energy is then obtained 
as an average of these two independent runs i^^'^^ 

In the case of the NVT ensemble the mass dependence 
of the Helmholtz free energy, _F, can be derived from the 
following relation analogous to Eq. ([7]) 



F{mp) = F{mH) 



{K{mi)) 



NVT 



dm J 



mi 



(11) 



D. Quantum-classical free energy difference 

The previous derivation can be applied to calculate the 
free energy difference of a classical system with respect to 
its quantum limit at the reference state point (Tr,Pr). 
The method implies the calculation of the reversible work 
associated to the following thermodynamic processes in 
the NPT ensemble: first (process A) , the mass of the wa- 
ter molecules in the quantum system is increased, so that 
one formally reaches the limit of infinite molecular mass, 
where the classical limit becomes exact; second (process 
B), the mass of the water molecules is reduced back to its 
actual value but now considering that the system remains 
in its classical limit. The reversible work performed along 
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these two processes is the Gibbs free energy difference be- 
tween the classical limit and the quantum system. This 
work can be calculated by Eq. ([7]) if applied to the partic- 
ular case where the masses of all the atoms in the system 
are scaled by the same factor. It it convenient to define 
the variable A as 



A 



Mo 
M 



(12) 



where Mq and M are the actual and the scaled molecular 
masses of water, respectively. The relative Gibbs free en- 
ergy difference, Gn^cia^ of the classical limit with respect 
to the quantum system is then 



phases over an overlapping range of temperatures. To 
this aim we consider the following thermodynamic rela- 
tion that relates the enthalpy, H{T), with the Gibbs free 
energy 



HiT) 



d [PG{T)] 
dp 



(15) 



where /? — l/fc^T is the inverse temperature. Here, we 
have omitted the explicit indication of the pressure and 
mass dependence of G and H , that are assumed to be 
constant. Integrating this equation between the inverse 
reference temperature, /3ij, and a final inverse tempera- 
ture, Pf, one gets 



GR,cla{Mo) 



\ \ A 



d\ . (13) 



The first summand in the integral determines the re- 
versible work of process A, while the second summand 
corresponds to process B. {K{M)) j^pj, is the ensemble 
average of the virial estimator for the total KE (sum 
of oxygen and hydrogen atoms) for the quantum system 
where the molecular mass of water is M = Mg/A. K^ia 
is the total KE in the classical limit that, as a conse- 
quence of the equipartition principle, is independent of 
the molecular mass M and depends only on tempera- 
ture {UbT /2 per degree of freedom). The integral in Eq. 

can be evaluated by the AS method, i.e., a PI sim- 
ulation is performed where the molecular mass is slowly 
increased as a function of A. The initial mass is Mq (set 
by A = 1) and the final mass is a large value defined by 
an adequate small Xp (typically around 0.01). The value 
of Garcia must be then obtained by an extrapolation of 
the integral to the value Xp 0, i.e.. 



G_R„ciQ = lim I{Xf) , 



(14) 



where I{Xf) is the definite integral in Eq. ([T^ when the 
upper integration limit is Xp instead of zero. In Appendix 
\X\ it is shown how the harmonic result for Eq. (|13p let us 
expect that the extrapolation can be accurately done by 
a simple polynomial fit in the variable Xp- An alterna- 
tive calculation of the free energy, Gji^cia, is the method 
of Morales and Singer,^'* based on a thermodynamic inte- 
gration that depends only on the potential energy. How- 
ever, this method is not so convenient for the calculation 
of free energy differences between isotopes, as it would 
require to use the classical limit as an intermediate step 
to set up a reversible path connecting both isotopes. 



ppGiTp) = I3bG{Tb 



dl3H{T) 



(16) 



/3r 



The RS method is implemented here by a non- 
equilibrium simulation where the inverse temperature is 
changed uniformly along a simulation run composed of J 
MDS as 



(17) 



where f5i is the inverse temperature of the i'th simu- 
lation step {i = 1,...,J) and the increment is A/3 = 
{f3p — /?fl)/( J — 1). The free energy at temperature Ti is 
discretized from Eq. (|16l) as 



AG(T,) = A-iG(7^,_i) + A/3 

'(18) 

H{Ti) is the system enthalpy at the i'th simulation step. 
This expression corresponds to the NPT ensemble and 
is also valid if G is substituted by the relative energy 
Gr, because f3G and I3Gr differ only in a temperature- 
independent constant, 5*0 [see Eq. ([1])]. The change of the 
inverse temperature Pi at each simulation step implies to 
update all variables that depend on temperature in our 
PI MD algorithm. These variables are the spring con- 
stants between beads, as well as the thermostat, barostat, 
and volume masses defined for the extended dynamics in 
the NPT ensemble.32-34 

A relation equivalent to Eq. (|18p can be derived in the 
NVT ensemble if G is substituted by the Helmholtz free 
energy, and the enthalpy, i?, by the internal energy 
U. 



III. TEST SIMULATIONS 



E. Temperature dependence of the free energy 



A. Radial distribution function 



The determination of the solid-liquid coexistence tem- 
perature requires to calculate the free energy of both 



The RDFs for 00, OH, and HH pairs have been cal- 
culated for water in the NVT ensemble at 298 K and 
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Figure 1: OO RDFs derived from quantum and classical sim- 
ulations of water at 298 K and density 0.997 g cm~^. For 
comparison the PI MD results of Ref. [12 are shown as open 
circles. 




Figure 3: HH RDFs derived from quantum and classical sim- 
ulations of water at 298 K and density 0.997 g cm^'^. For 
comparison the PI MD results of Ref. [l^ are shown as open 
circles. 




Figure 2: OH RDFs derived from quantum and classical sim- 
ulations of water at 298 K and density 0.997 g cm~^. For 
comparison the PI MD results of Ref. [12 are shown as open 
circles. 



Table I: Molecular properties (bond distance, bond angle and 
dipole moment) as well as kinetic (K) and potential energy 
(Upot) of liquid water at 298 K and density p = 0.997 g cm~^. 
ro - H is the distance of the goH RDF maximum associated 
to the H-bond. The KE is partitioned into H-isotope and O- 
atom contributions (Kj and Ko, respectively), ptmd is the 
maximum density of water at TMD, as derived from NPT 
simulations at P = 1 atm. Both classical and quantum results 
are given. The quantum results correspond to normal (H2O), 
heavy (D2O), and tritiated (T2O) water. 
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TMD (K) 


282(2) 






280(2) 
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1.004(2) 






1.002(2) 1.001(2) 



density 0.997 g cm~'^. The results of quantum and clas- 
sical simulations are presented in Figs. [T] to [31 The 
curves were derived from runs composed of 10^ MDS. 
Our quantum results have been compared to those pub- 
lished in Ref. il2j. We observe in the three RDFs an 
excellent agreement between both sets of independent 
calculations. The technical setup of the simulations dif- 
fers in many details, as the number of employed beads, 
the use of staging versus normal mode representations 
of the bead coordinates, the thermostats setup, the cut- 



off distance for short-range interactions, the approxima- 
tion used to perform the Ewald summation in reciprocal 
spacCfi^ and the schemes used for the RESPA molecular 
dynamics. Therefore, the agreement found by the RDFs 
curves provides a check for the accuracy of the employed 
computational conditions. 

The comparison of the RDFs obtained by PI simula- 
tions with the q-TIP4P/F model to experimental curves 
has been presented elsewhere^^ and will not be repeated 
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here. However, we focus here on differences found be- 
tween quantum and classical limits of the RDFs of this 
model. The 00 RDF curves in Fig. [l]show that classical 
water is more structured than in the quantum case. In 
the classical limit, the height of the first 00 peak is larger 
and the position of the maximum is displaced by about 
0.01 A towards shorter 00 distances. The second peak in 
the OH RDF curves in Fig. [2] corresponds to the H-bond 
distances. Controversial quantum and classical results 
have been published for this peak. A DFT study shows 
that in a quantum simulation this peak appears at shorter 
distances, which was interpreted as a hardening of the 
water structure with respect to classical simulations 4^ 
However, most simulations predict that quantum water is 
less structured than the classical counterpart.^ Our quan- 
tum results for the q-TIP4P/F model show that this peak 
appears at a distance 0.02 A larger than in the classical 
limit (see Table in agreement to the expectation that 
quantum corrections destabilize the H-bond network. 

A comparison of other equilibrium results of quantum 
and classical water simulations is summarized in Table 
m The quantum results have been derived for both nor- 
mal, heavy and tritiated water. We see that quantum 
corrections associated to the atomic mass increase the 
intramolecular OH bond length by more than 0.01 A, 
and the average molecular dipole moment, (/x), by 1.5%. 
Note that this increase in (/x) is expected to act against 
the destabilization of the H-bond network, as the electro- 
static interaction between neighboring water molecules 
becomes stronger. In Tab. |T] we also compare the av- 
erage kinetic and potential energies of water. We note 
that the KE of normal water is 5.45 kcal mol~^ larger 
in the quantum case. This increase is even larger for the 
potential energy (5.74 kcal mol~^) as a consequence of 
the anharmonicity of the model potential. A harmonic 
approximation predicts that both kinetic and potential 
energy increments must be identical. The partition of the 
KE between H- and 0-atom contributions shows that the 
three-fold rise of the KE of normal water in the quantum 
case is mainly due to the H-atoms. For comparison we 
note that ^°Ne is the Lennard- Jones-type system whose 
atomic mass is closest to that of a water molecule. For 
solid ^°Ne at 24 K (a temperature close to the melting 
point at atmospheric pressure) the rise in KE amounts to 
a factor of about 1.4 with respect to the classical limit. 



B. Temperature of maximum density 

We have performed quantum and classical NPT sim- 
ulations to calculate the temperature dependence of the 
water density, p, at ambient pressure. The results are 
summarized in Fig. Our quantum results, shown by 
open circles, are again in reasonable agreement to the PI 
MD data of Ref. [ij for the same model potential. The 
densities calculated by both sets of independent simula- 
tions differ by not more than 0.5%, which is of the order 
of the estimated statistical error. A systematic deviation 
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Figure 4: Density of water at 1 atm pressure obtained from 
classical and PI MD simulations. Lines are cubic polynomial 
fits to the simulation data. 



is found at temperatures corresponding to the maximum 
density, indicating that the convergence over density fluc- 
tuations is particularly difficult in this region. The TMD, 
as obtained by a cubic polynomial fit to the simulation 
points, is found at 280 K, in good agreement to previ- 
ously published results for this potential (see last two 
rows in Table |l|.^- A comparison of simulation results to 
experiment has been presented elsewhere?^ showing that 
the employed model provides realistic results. We recall 
that a proper description of the density of water is an im- 
portant requirement for a water model to account for hy- 
dration effects.*^ The classical TMD is found at slightly 
higher temperature than in the quantum case, the shift 
amounts to 2 K, which is of the order of the statisti- 
cal error. Classical and quantum temperature-density 
curves are so close that we have not tried to determine 
the hydrogen isotope effect in the TMD. However, we 
have checked that the density in simulations with tri- 
tiated water is, within the statistical error, very close 
to the normal water results. Therefore the employed q- 
TIP4P/F potential seems to be unable to reproduce the 
experimental isotope shift of the TMD in water. 

The small difference between the classical and quan- 
tum TMD is consistent with the previous work of Haber- 
shon et alJ^ who found similar trends by comparing clas- 
sical and quantum diffusion coefficients. A different re- 
sult was derived for a rigid water model, TIP4PQ/2005. 
For this model, the classical TMD was found about 27 K 
higher than the quantum temperature, and the isotope 
effect for tritiated water amounts to 16 K, overestimat- 
ing the experimental result of 9.4 K.^^ The main dif- 
ference between the rigid TIP4PQ/2005 and the flexible 
q-TIP4P/F potential is the presence of an anharmonic 
intramolecular potential in the latter. The mechanism 
that has been previously used to explain the striking 
differences between quantum and classical diffusion co- 
efficient of water when using either anharmonic or rigid 
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Table II: Relative free energies, Gr, of tritiated water at the 
reference point (Tr^Pu) as derived by independent AS simu- 
lations of different lengths. For a given simulation length two 
results are presented, corresponding to simulations where the 
initial and final integration limits {mH,mT) are interchanged. 
The last column shows the average of both independent runs. 



MDS 




Gr (kcal mol ^) 






rUH — >■ TTIT 


niT ^ rnn 


average 


10^ 


-5.947 


-5.950 


-5.949 


2x10^^ 


-5.951 


-5.949 


-5.950 


4x10^^ 


-5.949 


-5.951 


-5.950 


l.SxlO*' 


-5.950 


-5.952 


-5.951 



intramolecular modelsji^ can be also used here to explain 
the differences found in the quantum and classical TMDs 
when using either anharmonic or rigid models. This 
mechanism implies the coupling between anharmonic in- 
tramolecular stretches (OH) and intermolecular H-bonds 
(O- • • H) in water. Anharmonic zero point effects weaken 
both O- • • H and OH bonds (see the comparison of dis- 
tances {ron) and ro - H in Table HI quantum results are 
larger than classical ones). However, weakening of the 
OH bond increases the molecular dipole moment, that is 
accompanied by a concomitant strengthening of the in- 
termolecular H-bond. Therefore, when quantum effects 
associated to the nuclear masses are included in an an- 
harmonic flexible water potential, the H-bond is affected 
by competing factors acting in opposite directions, that 
might lead even to a mutual cancellation: a) weakening 
by zero point effects; b) strengthening by the increase of 
the molecular dipole moment. Obviously this competing 
mechanism is absent if the water model is rigid. Which 
of the two competing effects is dominant will depend on 
the details of the potential model and on the physical 
property under consideration. For the TMD of the q- 
TIP4P/F model, it seems that both competing effects 
cancel each other, leading to nearly identical quantum 
and classical curves of the liquid density as a function 
of temperature and then to a vanishingly small isotope 
effect in the TMD of water. 



IV. ISOTOPE EFFECTS IN THE MELTING 
TEMPERATURE 

A. D2O and T2O 

To study the isotope effect in the melting temperature 
of ice we need first to calculate the relative free energies, 
Gr, of deuterated and tritiated phases at the reference 
state point (Tr^Ph). To this aim we have performed 
AS simulations in the NPT ensemble according to Eq. 
(fTO| . For the liquid phase, the dependence of Gr with 
the isotope mass is shown in Fig. [Sj The convergence 
of Gr has been checked by comparing independent non- 
equilibrium simulations of different lengths for both solid 
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Figure 5: Relative free energy of liquid water as a function 
of the hydrogen isotope mass. The open circles corresponds 
to the masses of H, D, and T, respectively. The result was 
obtained by non-equilibrium simulations with the AS method 
at the reference state point {Tr,Pr). 



Table III: Relative free energy, Gr, at the reference point 
{Tr = 251 K, Pr= 1 atm) of solid and liquid phases of normal, 
heavy, and tritiated water. Gr is given in kcal mol~^, and 
its estimated error is ±0.001 kcal mol~^. The last column 
summarizes the results obtained in the classical limit {Gr^cIo)- 





H2O 


D2O 


T2O 


classical limit 


solid (s) 





-4.108 


-5.987 


-8.912 


liquid (1) 





-4.079 


-5.951 


-8.875 


s - 1 





-0.028 


-0.036 


-0.036 



and liquid phases. This test is presented in Table Ull for 
tritiated water. We find that the statistical error in the 
free energy is even lower for heavy water and for the solid 
phases. Even with a modest number of simulation steps 
(J = 10^) the AS method shows reasonable convergence. 
Our final results for Gr are summarized in Table iHll cor- 
responding to non-equilibrium simulations with 1.5x10^ 
MDS for the liquid and 4 x 10^ MDS for the sohd phases. 
For the heavier isotopes (D2O and T2O) the relative free 
energy of the solid at the reference state point is lower 
than in the liquid phase. This difference determines an 
isotope shift in the melting temperature, as temperature 
must be increased to restore the coexistence condition of 
equal free energy for both phases. The free energy differ- 
ence between solid and liquid is larger for tritium than 
for deuterium, therefore the isotope shift in the melting 
temperature is expected to be larger for tritiated water 
than for deuterated water. 

The temperature dependence of Gr_ has been cal- 
culated by RS simulations in the NPT ensemble [see 
Eq. p^ ]. Several independent simulations of different 
lengths were performed to fix the conditions of adequate 
convergence for the non-equilibrium simulations. Again, 
each free energy estimation was based on two indepen- 
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Figure 6: Relative free energy of normal water and ice at 
pressure of 1 atm as determined by our RS simulations. The 
melting point is Tm- 



Figure 8: Relative free energy of tritiated water and ice at 
pressure of 1 atm as determined by our RS simulations. The 
melting point is Tm- 
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Table IV: Computational conditions used in the non- 
equilibrium RS simulations to determine the temperature de- 
pendence of the relative Gibbs free energy of the studied iso- 
tope compositions of water and ice. 



isotope 


phase 


T range (K) 


MDS 


H2O 


liquid 


250-252 


2x10^ 


H2O 


solid 


250-252 


5x10^* 


D2O 


liquid 


251-263 


2x10** 


D2O 


solid 


251-259 


2x10^5 


T2O 


liquid 


251-263 


2x10** 


T2O 
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solid 
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Figure 7: Relative free energy of deuterated water and ice at 
pressure of 1 atm as determined by our RS simulations. The 
melting point is Tm- 



dent simulations where the integration limits of the tem- 
perature in the reversible path were interchanged, and 
the final results were obtained as an average of both 
RS simulations. The employed computational conditions 
(integration limits for the temperature and number of 
MDS) used in our RS simulations are summarized in 
Table |IVl The liquid RS simulations require typically 
simulation lengths one order of magnitude larger than 
the solid phase. The Gii{T) results obtained for normal, 
deuterated and tritiated phases are presented in Figs. [6]- 
ISl The condition of equal free energy of the solid and 
liquid phases determines the melting point at ambient 
pressure. The melting temperature is shifted towards 



higher values as the hydrogen isotope mass increases. 
The estimated isotope shift in the melting temperature 
is 6.5±0.5 K for heavy water and 8.2±0.5 K for tritiated 
water for the q-TIP4P/F model. Experimental isotope 
shifts amount to 3.8 K and 4.5 K, respectively. Thus, the 
q-TIP4P/F model overestimates the isotope shift in the 
melting temperature of ice. This result contrasts with 
the nearly absence of isotope shifts in the TMD of water, 
as predicted by the same potential. 

Melting in atomic solids can be thought of as being 
controlled by the local vibrations of the atoms. The Lin- 
demann criterion presents a threshold for the maximum 
amplitude of atomic vibrations that can be sustained by 
the crystal. At the melting point of ice we find that 
this amplitude is similar for both O- and H-atoms and 
amounts to about 6% of the H-bond distance. 
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Figure 9: The function I{Xp) obtained for liquid water by 
non-equilibrium AS NPT simulations at the reference state 
point {ThiPr)- The extrapolation Af — 5- is shown for the 
solid and liquid phases in the inset of the figure. 
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Figure 10: Relative free energy of water and ice at pressure 
of 1 atm as determined by our RS simulations in the classical 
limit. The melting point is Tm- 



B. Classical limit 

The determination of the classical melting point im- 
plies to calculate the relative free energy, Garcia-, of the 
classical limit of water and ice with respect to the quan- 
tum case at the reference state point. To this aim we have 
performed AS simulations where the molecular mass of 
water is slowly changed from its normal value, Mq, to an 
arbitrary large value, Mp, determined by the parameter 



Xp = Mo/Mp. The function I{Xp) introduced in Sec. 
IIIDI [see Eq. is then obtained by a discretization 

of the integral given in Eq. (|13p as a running average 
along the non-equilibrium simulation. The AS results 
in the solid phase have been derived with 4x10^^ MDS 
while the liquid phase required longer runs of 1.6x10^ 
MDS. As usual in our non-equilibrium simulations, two 
independent simulation runs were performed by inter- 
changing the integration limits, and the final free energy 
is the average of these two runs. The function I{Xf) for 
the liquid phase is shown in Fig. [HI The inset in the 
figure shows the result of the extrapolation Xp — > for 
both solid and liquid phases. The extrapolated values 
are summarized in the last column of Table IIIIl These 
figures were obtained by a polynomial fit of /(Af), af- 
ter checking that the extrapolation is reasonably stable 
against both the polynomial degree and the Xp interval 
employed in the numerical fit. The results given in Table 
mil were derived by a lO'th degree polynomial fit in the 
Xp range [0.01,0.2]. Somewhat surprising is that the dif- 
ference between the relative free energy of the solid and 
liquid phases in the classical limit is, within the statisti- 
cal error, identical to the value found for tritiated phases. 
This fact will be relevant for the melting temperature of 
the classical system. 

The temperature dependence of the relative free energy 
of water and ice at ambient pressure is presented in Fig. 
[TO] for the classical limit. The computational conditions 
for the non-equilibrium RS simulations were summarized 
in the last two rows of Table IIVI The estimated classical 
melting point is 259.3±0.5 K. We find again an excellent 
agreement to the result reported by Habershon et alJ^ of 
259±1K. Note that both classical melting temperatures 
were determined by different methods, i.e., by direct co- 
existence simulationsi^ and by free energy calculations. 
An advantage of the free energy method is that other im- 
portant physical quantities are readily available, which is 
not the case by a direct coexistence method. 

In Table |V] we summarize the simulation results of 
several melting properties for the studied isotopic com- 
positions of water and also in the classical limit. The 
melting entropy, ASm, was estimated from the numeri- 
cal temperature derivative of the Gibbs free energy curves 
given in Figs. [51151 and [TOl at coexistence conditions, while 
the melting enthalpy was calculated by two independent 
ways: (a) as Tm^Sm', (b) by direct calculation of the 
solid and liquid enthalpies by NPT simulations at coex- 
istence conditions. The agreement between both meth- 
ods provides evidence about the internal consistency of 
our results. The isotope effect in the melting entropy 
and enthalpy is low, practically within the statistical un- 
certainty of our results. The employed water model un- 
derestimates both the experimental values of the melting 
entropy and enthalpy. At coexistence conditions the neg- 
ative sign in the change of the KE upon melting indicates 
that the KE of ice is larger than that of the liquid phase. 
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Table V: Melting temperature, entropy and enthalpy for normal, heavy, and tritiated water as well as classical limit results at 
ambient pressure. The melting enthalpy was estimated by two independent methods. The change in the kinetic and potential 
energy upon melting (liquid minus solid values) and the molar volume (Kn) of the solid and liquid phases are also given. The 
standard error in the final digits is given in parenthesis. 





classical 


T2O 


D2O 


H2O 


exp. H2O 


Tm{K) 


259.3(5) 


259.2(5) 


257.5(5) 


251 


273.15 


AS™ (cal mor^K^^) 


4.52(5) 


4.45(5) 


4.52(5) 


4.40(5) 


5.3 


AHm" (kcal mor^) 


1.17(2) 


1.15(2) 


1.16(2) 


1.10(2) 


1.44 


AHm'' (kcal mor^) 


1.18(3) 


1.10(3) 


1.13(3) 


1.07(3) 


1.44 


AKm ( kcal mor^ ) 





-0.024(5) 


-0.038(5) 


-0.078(5) 




AL'pot(kcal moP^) 


1.18(3) 


1.12(3) 


1.17(3) 


1.14(3) 




Kn,s(cm^mol~^) 


19.40(1) 


19.47(1) 


19.49(1) 


19.56(1) 


19.66 


Kn,i(cm^mor^) 


17.99(6) 


18.15(6) 


18.03(6) 


17.96(6) 


18.02 


AKr,(cm^mor^) 


-1.41(6) 


-1.32(6) 


-1.46(6) 


-1.60(6) 


-1.64 



"from TmAS m 

'from independent solid and liquid NPT simulations at coexis- 
tence conditions 



At room temperature the KE of tritium in tritiated wa- 
ter is a factor about 2.4 times larger than in the classical 
limit (see Table i.e., although tritium is three times 
heavier than hydrogen, quantum effects related to its nu- 
clear mass are still significant below room temperature. 
Therefore an unexpected result of the employed model 
potential is that the classical melting point is nearly iden- 
tical to that one found for tritiated water, and the reason 
for this behavior will be investigated in the next subsec- 
tion. 



C. Kinetic energy and molecular mass 

We have already seen that isotope shifts in the melting 
point are caused by the fact that the free energy of the 
solid and liquid phases changes by different amounts as 
the isotopic mass changes. This is a quantum effect re- 
lated to the atomic mass, as in the classical limit the mass 
dependence of the free energy is identical for both phases. 
To understand the shift found in the melting temperature 
in the classical limit, it is interesting to study the differ- 
ence between the Gibbs free energy of the solid and liquid 
phases as a function of the scaled molecular mass, Mp, 
as determined by the parameter Xp = Mo/Mp. This 
difference is given by an expression similar to Eq. (jl3p 



GrAMp) 



Gr^Mp) = 
(jr.(M))^p^ 
A 



{Ki{M)) 



NPT 



A 



dX . (19) 



The instantaneous values of the virial KE estimators of 
the integrand are available from the non-equilibrium AS 
simulations used to calculate Fig. |9l These values have 
been used to plot the solid-liquid Gibbs free energy dif- 
ference as a function of the inverse molecular mass, Mp^ , 
in Fig. 111! This curve is not a monotonous function but 
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Figure 11: Gibbs free energy difference between ice and water 
as a function of the molecular mass at the reference state 
point {Tr,Pr). Depending on the molecular mass the KE of 
the liquid {Ki) may be larger than that of the solid {Ks). 



it displays a minimum for a molecular mass defined by 
the parameter Xp ^ 0.14. This minimum makes that the 
free energy difference in the classical limit (corresponding 
to Xp — 0, open circle in Fig. [TT|) is similar to that ob- 
tained for tritiated phases (see Table IIIip , and therefore 
the classical melting point is found at nearly the same 
temperature as for T2O. If the curve in Fig. [Tl]were a 
monotonous decreasing function up to Xp = (i.e., with- 
out the presence of a minimum) , then the classical abso- 
lute value of the free energy difference between the two 
phases should be larger than for T2O, and therefore the 
classical melting temperature should also increase with 
respect to the tritiated phase. 
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Figure 12: Kinetic energy difference between ice and water as 
a function of the molecular mass at the reference state point 
(Th^Ph). The line is a guide to the eye. 



The existence of a minimum in the function given in 
Fig. [11] is related to a change in the sign of the integrand 
of Eq. p9p in the region around the minimum. This fact 
is clearly seen by plotting in Fig. [T^] the difference be- 
tween the KE of solid and liquid phases as a function of 
the inverse molecular mass, Mp^, at the reference state 
point. We find that when the molecular mass is large 
(Af < 0.14, or > 7Mo ) the KE of water is larger 
than that of ice, while the opposite behavior is found for 
lower molecular masses {Xp > 0.14 ). The physical ori- 
gin of this appareirtly complicate behavior of the KE is 
related to the presence of the two types of bonds in the 
water phases: intramolecular OH pairs and intermolecu- 
lar H-bonds. The H-bond network is stronger in ice than 
in water, as a result of the higher molecular disorder in 
the liquid, which implies that vibrational frequencies of 
H-bonds are larger in ice than in water. On the con- 
trary, the OH stretch frequency in the liquid is higher 
than in the solid. In Fig. [T3]we display the OH distance 
for water and ice as a function of the molecular mass at 
the reference state point. For all molecular masses, the 
OH distance is lower for water than for ice, a fact that 
implies also a higher vibrational OH stretch frequency in 
water. This behavior is illustrated by the inset in Fig. [T3] 
that displays the increase of the quasi-harmonic stretch 
frequency ujqh for the employed q-TIP4P /F potential as 
the OH distance decreases. 

In the classical limit, corresponding to Af = (infinite 
molecular mass), all modes have the same KE (equipar- 
tition principle) and the KE difference between ice and 
water vanishes. In the case of large, but finite, molecular 
mass {0 < Xp < 0.14) , the KE differeirce is determiired 
by the modes with highest vibrational frequencies, i.e., 
the OH stretches. The reason is that the leading quan- 
tum correction for vibrational modes, as the molecular 
mass decreases, corresponds to those modes that satisfy 
that their energy quantum is larger than the thermal en- 
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Figure 13: Intramolecular OH distance of ice and water as 
a function of the molecular mass at the reference state point 
(Tr^Pr). Lines are guides to the eye. The inset shows the 
quasi-harmonic stretch frequency looh (in cm~^) as a function 
of the OH distance for the q-TlP4P/F potential, uuoh was 
derived from the second derivative of the potential energy 
with respect the OH bond distance and by considering the 
actual O and H masses. 



ergy, fiw > fc^Tfl, and this condition will be first met for 
the modes with highest frequency (OH stretches). Thus, 
for large molecular masses with Xp < 0.14, the KE of wa- 
ter results larger than that of ice, as the stretch frequen- 
cies ujoH are larger in the liquid. For smaller molecular 
masses {Xp > 0.14), the quantum behavior of the vibra- 
tional modes associated to H-bonds becomes important. 
For H-bonds related modes, the KE is larger in ice than 
in water, as their vibrational frequencies are higher in 
ice. The effect of the H-bonds dominates over the OH 
stretches for Xp > 0.14 in the sense that the KE of ice 
results larger than that of water. 

The shorter OH distance in liquid water is a conse- 
quence of the coupling between OH stretches and H- 
bonds. Since the H-bond network in water is weaker 
than in ice, then the intramolecular OH bonds become 
shorter. A recent Compton scattering study of water ver- 
sus ice Ih arrives to the same conclusion: an elongation of 
the H-bond in water leads to a systematic shortening of 
the intramolecular OH bond, and the OH bond in water 
is about 0.01 A shorter than in ice4^ 



V. CONCLUSIONS 

In the present work we have used free energy tech- 
niques to study the isotope shift in the meltiirg temper- 
ature of ice Ih. The employed AS and RS methods are 
based on algorithms where the Hamiltonian or a state 
variable are adiabatically changed along a simulation run. 
The reversible work associated to this change is equal to 
the free energy difference between the initial and final 
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state, as long as the change is performed slowly (adia- 
batically) . These free energy algorithms are easily imple- 
mented in any code prepared for equilibrium simulations. 

Our PI simulations of water have been done with the q- 
TIP4P/F model, a point charge potential that includes 
an anharmonic treatment of molecular flexibility. Sev- 
eral equilibrium properties (RDFs, liquid density as a 
function of temperature) have been compared to results 
previously published for this potential. The agreement 
found between independent simulations provides a check 
for the employed computational method and simulation 
conditions. 

We have found that the experimental isotope effect in 
the TMD of water is not correctly described by the em- 
ployed potential. In fact, the TMD calculated in the 
classical limit is very close to the quantum value. Quan- 
tum effects associated to the atomic mass activate a com- 
peting mechanism that produces a weakening of the H- 
bond through quantum zero point fluctuations, but also 
a strengthening of the same H-bonds by the increase in 
the molecular dipole moment of water molecules, that is 
found in the quantum simulations. Although this com- 
peting mechanism, previously observed by Habershon et 
al}^ to explain the trends in the classical and quantum 
diffusion coefHcients, seems to be a real physical effect, 
its influence in a given physical property depends on the 
details of the anharmonic flexible water potential, so that 
a small imbalance between both competing factors might 
lead to an unphysical result. This seems to be the case 
for the isotope effect in the TMD of water. 

The isotope effect in the melting temperature at am- 
bient pressure of ice has been determined by the calcula- 
tion of the Gibbs free energy as a function of the isotope 
mass and temperature. We find that the isotope effect 
predicted by the q-TIP4P/F potential (6.5±0.5 K and 
8.2±0.5 K for heavy and tritiated water, respectively) 
are larger than the experimental values (3.8 K and 4.5 
K, respectively). An unexpected result is that the classi- 
cal melting point at 1 atm is nearly identical to the one 
obtained for tritiated water. We have shown that this 
behavior is related to the fact that the OH stretches in 
water display a higher frequency than in the solid phase. 
It is clear that the employed q-TIP4P /F potential is not 
able to quantitatively predict the isotope effect in the 
melting temperature of ice, but it has helped to identify 
that the coupling between molecular fiexibility and the 
H-bond network may have significant implications in the 
phase behavior of water. 
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Appendix A: Extrapolation of Gn^cia 

The largest deviation between a quantum and classi- 
cal treatment of a water phase (either solid or liquid) is 
expected to be due to those vibrational modes of high- 
est frequency (i.e., intramolecular OH stretches and H- 
bonds) . A harmonic treatment of these modes gives use- 
ful information on the analytical behavior of the definite 
integral in Eq. (jl9p when the upper integration limit 
vanishes (Af — >■ 0). Note that the limit Xp = is not ac- 
cessible numerically as it would correspond to an infinite 
mass. Let us assume a simple one-dimensional harmonic 
oscillator defined by a wavenumber wq and mass Mq. The 
oscillator mass depends on the integration parameter A 
as M = Mq/A, while the oscillator wavenumber varies as 
Lo = wqA"'^/^. The thermal expectation value of the KE of 
the quantum oscillator for a given A is given by 



{K)=K,\^/^coth(^\^l-A , (Al) 



where K^ia is the classical KE 



Krin. ^ -knT 



and Kq is the zero-point KE 



Ko = — 



(A2) 



(A3) 



By a Taylor expansion of the r.h.s. of Eg. ljAip . the inte- 
grand of Eq. (jl3p for a harmonic mode can be written 
as 



(K) - K,la 

A 



Kcla 



^cla 



A-f 0(A2) 



(A4) 



The definite integral given in Eq. must be obtained 
by numerical extrapolation of I{Xf) to the Xp — limit 
[see Eq. p^ ]. The harmonic result for the integrand in 
Eq. (jA4p suggests that a polynomial fit for I{Xf) is a 
convenient numerical way to obtain the desired extrapo- 
lation. As the employed q-TIP4F/P potential is anhar- 
monic, we do not expect to assign any particular meaning 
to the fitted coefficients. 
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